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We present a perturbative method for calculating phonon properties of an insulator in the pres- 
ence of a finite electric field. The starting point is a variational total-energy functional with a 
field-coupling term that represents the effect of the electric field. This total-energy functional is 
expanded in small atomic displacements within the framework of density-functional perturbation 
theory. The linear response of field-polarized Bloch functions to atomic displacements is obtained 
by minimizing the second-order derivatives of the total-energy functional. In the general case of 
nonzero phonon wavevector, there is a subtle interplay between the couplings between neighboring 
k-points introduced by the presence of the electric field in the reference state, and further-neighbor 
k-point couplings determined by the wavevector of the phonon perturbation. As a result, terms 
arise in the perturbation expansion that take the form of four-sided loops in k-space. We implement 
the method in the ABINIT code and perform illustrative calculations of the field-dependent phonon 
frequencies for III-V semiconductors. 

PACS numbers: 63.20.Dj, 78.20.Jq, 63.20.Kr, 71.55.Eq, 71.15.-m, 77.65.-j. 
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I. INTRODUCTION 

The understanding of ferroelectric and piezoelectric 
materials, whose physics is dominated by soft phonon 
modes, has benefited greatly from the availability of first- 
principles methods for calculating phonon properties. In 
general, these methods can be classified into two main 
types, the direct or frozen-phonon approach 1 ^ and the 
linear-response approach^^ In the former approach, the 
properties of phonons at commensurate wavevectors are 
obtained from supercell calculations of forces or total- 
energy changes between between equilibrium and dis- 
torted structures. In the latter approach, based on 
density- functional perturbation theory (DFPT), expres- 
sions are derived for the second derivatives of the total 
energy with respect to atomic displacements, and these 
are calculated by solving a Sternheimer equation^ or by 
using minimization methods. 4 ' 5 Compared to the direct 
approach, the linear-response approach has important 
advantages in that time-consuming supercell calculations 
are avoided and phonons of arbitrary wavevector can be 
treated with a cost that is independent of wavevector. 
However, existing linear-response methods work only at 
zero electric field. 

The development of first-principles methods for treat- 
ing the effect of an electric field 8 in a periodic system has 
been impeded by the presence of the electrostatic poten- 
tial £ r in the Hamiltonian. This potential is linear in real 
space and unbounded from below, and thus is incompat- 
ible with periodic boundary conditions. The electronic 
bandstructure becomes ill-defined after application of a 
potential of this kind. Many attempts have been made 
to overcome this difficulty. For example, linear-response 
approaches have been used to treat the electric field as 
a perturbationi 5 ^ It is possible to formulate these ap- 
proaches so that only the off-diagonal elements of the 
position operator, which remain well defined, are needed, 



thus allowing for the calculation of Born effective charges, 
dielectric constants, etc. Since it is a perturbative ap- 
proach, a finite electric field cannot be introduced. 

Recently, a total-energy method for treating insula- 
tors in nonzero electric fields has been proposedi^ In 
this approach, an electric enthalpy functional is defined 
as a sum of the usual Kohn-Sham energy and an £ ■ P 
term expressing the linear coupling of the electric field to 
the polarization P. The enthalpy functional is minimized 
with respect to field-polarized Bloch states, and the infor- 
mation on the response to the electric field is contained 
in these optimized Bloch states. Using this approach, it 
is possible to carry out calculations of dynamical effec- 
tive charges, dielectric susceptibilities, piezoelectric con- 
stants, etc., using finite-difference methodsi*^ It would 
also be possible to use it to study phonon properties in 
finite electric field, but with the aforementioned limita- 
tions (large supercells, commensurate wavevectors) of the 
direct approach. 

In this work, we build upon these recent developments 
by showing how to extend the linear-response methods 
so that they can be applied to the finite-field case. That 
is, we formulate DFPT for the case in which the unper- 
turbed system is an insulator in a finite electric field. 
Focusing on the case of phonon perturbations, we de- 
rive a tractable computational scheme and demonstrate 
its effectiveness by carrying out calculations of phonon 
properties of polar semiconductors in finite electric fields. 

This paper is organized as follows. In Sec. [H] we review 
the total-energy functional appropriate for describing an 
insulator in an electric field, and discuss the effect of the 
electric field on the phonon frequencies both for our ex- 
act theory and for a previous approximate theory. The 
second-order expansion of the total-energy functional is 
derived in Sec. IIIII and expressions for the force-constant 
matrix are given, first for phonons at the Brillouin zone 
center and then for arbitrary phonons. In Sec. lIVI we re- 
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port some test calculations of field-induced changes of 
phonon frequencies in the III-V semiconductors GaAs 
and AlAs. Section V contains a brief summary and con- 
clusion. 



II. BACKGROUND AND DEFINITIONS 

A. Electrical enthalpy functional 

We start from the electric enthalpy functional 7 

F[R; ft £] = E KS [R; ft-Q£- P mac [ft, (1) 

where -Eks has the same form as the usual Kohn-Sham 
energy functional in the absence of an electric field. Here 

is the cell volume, P mac is the macroscopic polariza- 
tion, £ is the homogeneous electric field, R are the atomic 
positions, and ip are the field-polarized Bloch functions. 
Note that P mac has both ionic and electronic contribu- 
tions. The former is an explicit function of R, while 
the latter is an implicit function of R through the Bloch 
functions, which also depend on the atomic positions. 
When an electric field is present, a local minimum of this 
functional describes a long-lived metastable state of the 
system rather than a true ground state (indeed, a true 
ground state does not exist in finite electric field) i 

According to the modern theory of polarization, 9 the 
electronic contribution to the macroscopic polarization is 
given by 

ief M f 

Pmac = , To ^ / dk(M kn |Vfe|M k „), (2) 
W n=l j BZ 

where e is the charge of an electron (e < 0), f=2 for 
spin degeneracy, M is the number of occupied bands, Mkn 
are the cell-periodic Bloch functions, and the integral is 
over the Brillouin zone (BZ). Making the transition to a 
discretized k-point mesh, this can be written in a form 




that is amenable to practical calculations. In this expres- 
sion, for each lattice direction i associated with primitive 
lattice vector a^, the BZ is sampled by N± strings of k- 
points, each with Ni points spanning along the reciprocal 
lattice vector conjugate to a^, and 

(Skk')mn = (u m k\Unk') (4) 

are the overlap matrices between cell-periodic Bloch vec- 
tors at neighboring locations along the string. Because 
Eqs. (12i;St can be expressed in terms of Berry phases, this 
is sometimes referred to as the "Berry-phase theory" of 
polarization. 



B. Effect of electric field on phonon frequencies 

1. Exact theory 

We work in the framework of a classical zero-temper- 
ature theory of lattice dynamics, so that quantum zero- 
point and thermal anharmonic effects are neglected. In 
this context, the phonon frequencies of a crystalline insu- 
lator depend upon an applied electric field in three ways: 
(i) via the variation of the equilibrium lattice vectors (i.e., 
strain) with applied field; (ii) via the changes in the equi- 
librium atomic coordinates, even at fixed strain; and (iii) 
via the changes in the electronic wavefunctions, even at 
fixed atomic coordinates and strain. Effects of type (i) 
(essentially, piezoelectric and clcctrostrictivc effects) are 
beyond the scope of the present work, but are relatively 
easy to include if needed. This can be done by computing 
the relaxed strain state as a function of electric field using 
the approach of Ref. and then computing the phonon 
frequencies in finite electric field for these relaxed struc- 
tures using the methods given here. Therefore, in the 
remainder of the paper, the lattice vectors are assumed 
to be independent of electric field unless otherwise stated, 
and we will focus on effects of type (ii) ( "lattice effects" ) 
and type (iii) ( "electronic effects" ) . 

In order to separate these two types of effects, we first 
write the change in phonon frequency resulting from the 
application of the electric field as 

A0(q;£)=0(q;R £ ,£)-0(q;R o ,O), (5) 

where 0(q; R, £ ) is the phonon frequency extracted from 
the second derivative of the total energy of Eq. JIJ with 
respect to the phonon amplitude of the mode of wavevec- 
tor q, evaluated at displaced coordinated R and with 
electrons experiencing electric field £. Also, Rf are the 
relaxed atomic coordinates at electric field £ , while Ro 
are the relaxed atomic coordinates at zero electric field. 
Then Eq. (J5J can be decomposed as 

A0(q;£) = A0 cl (q;£)+A0 ion (q;£) (6) 

where the electronic part of the response is defined to be 

A 0cl (q; £) = 0(q; R , £) - 0(q; R , 0) (7) 

and the lattice (or "ionic" ) part of the response is defined 
to be 

A 0ion (q; £) = 0(q; R £ , £ ) - (q; R , £). (8) 

In other words, the electronic contribution reflects the 
influence of the electric field on the wavefunctions, and 
thereby on the force-constant matrix, but evaluated at 
the zero-field equilibrium coordinates. By contrast, the 
ionic contribution reflects the additional frequency shift 
that results from the field- induced ionic displacements. 

The finite-electric-field approach of Refs. 00 provides 
the methodology needed to compute the relaxed coordi- 
nates Rg , and the electronic states, at finite electric field 
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£. The remainder of this work is devoted to developing 
and testing the techniques for computing 0(q;R, £) for 
given q, R, and £, needed for the evaluation of Eq. 
We shall also use these methods to calculate the vari- 
ous quantities needed to perform the decomposition of 
Eqs. (ItiliSll . so that we can also present results for A0 e i 
and A0ion separately in Sec. IIVI 



2. Approximate theory 

Our approach above is essentially an exact one, in 
which Eq. (0 is evaluated by computing all needed quan- 
tities at finite electric field. However, we will also com- 
pare our approach with an approximate scheme that 
has been developed in the literature over the last few 
years, 10 ' 11-12 ' 13 in which the electronic contribution is 
neglected and the lattice contribution is approximated 
in such a way that the finite-electric-field approach of 
Refs. 0H is not needed. 

This approximate theory can be formulated by starting 
with the approximate electric enthalpy functional^ 



F[R; £] 



<[R]-^-P<l[R], 



(9) 



where [R] is the zero-field ground-state Kohn-Sham 

energy at coordinates R, and PmL is the corresponding 
zero-field electronic polarization. In the presence of an 
applied electric field £ , the equilibrium coordinates that 
minimize Eq. satisfy the force-balance equation 



(0) 
KS 



rfR 



Z<°> • £ = 



(10) 



where = fl c?Pmac/dR is the zero-field dynamical ef- 
fective charge tensor. That is, the sole effect of the elec- 
tric field is to make an extra contribution to the atomic 
forces that determine the relaxed displacements; the elec- 
trons themselves do not "feel" the electric field except in- 
directly through these displacements. In Ref. 10, it was 
shown that this theory amounts to treating the coupling 
of the electric field to the electronic degrees of freedom 
in linear order only, while treating the coupling to the 
lattice degrees of freedom to all orders. Such a theory 
has been shown to give good accuracy in cases where the 
polarization is dominated by soft polar phonon modes, 
but not in systems in which the electronic and lattice 
polarizations are comparable>i2iii*i2iiii4 

In this approximate theory, the effect of the electric 
field on the lattice dielectric properties^ and phonon 
frequencies^ comes about through the field-induced 
atomic displacements. Thus, in the notation of Eqs. 
IHJ), the frequency shift (relative to zero field) is 

A ; on (q; £) = 0(q; R^, 0) - 0(q; R , 0) (11) 

in this approximation, where R^ is the equilibrium po- 
sition according to Eq. I|10|) . We will make compar- 
isons between the exact R^ and the approximate Rg, 



and the corresponding frequency shifts A0i on (q, £) and 
A0| on (q,£) later in Sec. IIVI 



III. PERTURBATION EXPANSION OF THE 
ELECTRIC ENTHALPY FUNCTIONAL 



We consider an expansion of the properties of the sys- 
tem in terms of small displacements A of the atoms away 
from their equilibrium positions, resulting in changes in 
the charge density, wavefunctions, total energy, etc. We 
will be more precise about the definition of A shortly. We 
adopt a notation in which the perturbed physical quan- 
tities are expanded in powers of A as 



Q(A) 



)(0) 



+ AQ« + A 2 Q< 2) + A 3 



)(3) 



+ .. 



(12) 



where Q( n ' = [\/n\)d n Q / d\ n . The immediate depen- 
dence upon atomic coordinates is through the external 
potential v ox t(A), which has no electric- field dependence 
and thus depends upon coordinates and pseudopotentials 
in the same way as in the zero-field case. The changes 
in electronic wavefunctions, charge density, etc. can then 
be regarded as being induced by the changes in w C xt- 



A. Zero q wavevector case 



The nuclear positions can be expressed as 



R„ 



(13) 



where t„ is a lattice vector, d„ is a basis vector within 
the unit cell, and b„„ is the instantaneous displacement 
of atom v in cell n. We consider in this section a phonon 
of wavevector q = 0, so that the perturbation does not 
change the periodicity of the crystal, and the perturbed 
wavefunctions satisfy the same periodic boundary condi- 
tion as the unperturbed ones. To be more precise, we 
choose one sublattice v and one Cartesian direction a 
and let b nva = A (independent of rt), so that we are 
effectively moving one sublattice in one direction while 
while freezing all other sublattice displacements. Since 
the electric enthalpy functional of Eq. Q is variational 
with respect to the field-polarized Bloch functions under 
the constraints of orthonormality, a constrained varia- 
tional principle exists for the second-order derivative of 
this functional with respect to atomic displacements. 15 
In particular, the correct first-order perturbed wavefunc- 
tions V'mk can be- obtained by minimizing the second- 
order expansion of the total energy with respect to atomic 
displacements, 



(2) 



-fiP ma c[^;^]-£) , (14) 
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subject to the constraints 



(15) 



(where m and n run over occupied states). The fact that 
only zero-order and first-order wavefunctions appear in 
Eq. Ijl4(l is a consequence of the "2n+l theorem."^ 

Recalling that IV'ik) ^ s * ne first-order wavefunction re- 
sponse to a small real displacement A of basis atom v 
along Cartesian direction a, we can expand the external 
potential as 



w«t(r) = v£>(r) + ^xLa(r)A + w£L a (r)A a 



(16) 



terms, respectively. The first and last of these are given 
by 

j, occ 

FkS = -TJ ^2(lpnk\T + WextlV'nk) + ^HxcM, (21) 



kn 



and 



j, occ 

-FLM = —Jj A ki „ m ((V> mk |V'nk) - <?>mn), (22) 

k,mn 

where iV is the number of k-points in the BZ. As for 
the Berry-phase term, we modify the notation of Eq. (|3J) 
slightly to write this as 



where 



,(!) 



E 



dfcxt(r) 



(17) 



where 



(23) 



,(2) 



,» = E 



dVxt(r) 
dR 2 



(18) 



etc. From this we shall construct the second-order energy 
F^ 2 ' of Eq. 114(1. which has to be minimized in order 

to find IV'ik)- The minimized value of F^> gives, as 
a byproduct, the diagonal element of the force-constant 
matrix associated with displacement va. Once the IV'ik) 
have been computed for all va, the off-diagonal elements 
of the force-constant matrix can be calculated using a 
version of the 2n + 1 theorem as will be described in 
Sec. 1111 A 31 



1. Discretized k mesh 

In practice, we always work on a discretized mesh of k- 
points, and we have to take into account the orthogonal- 
ity constraints among wavefunctions at a given k-point 
on the mesh. Here, we are following the "perturbation 
expansion after discretization" (PEAD) approach intro- 
duced in Ref. Il6l That is, we write down the energy 
functional in its discretized form, and then consistently 
derive perturbation theory from this energy functional. 
Introducing Lagrange multipliers A krnn to enforce the 
orthonormality constraints 



(VVriklVVik) = S mn , 



(19) 



where VVik are the Bloch wavefunctions, and letting N be 
the number of k-points, the effective total-energy func- 
tional of Eq. can be written as 



F = F, 



KS 



Fi 



BP 



LM 



(20) 



where F K s = E KS , F B p 



-HP„ 



■ £, and Flm are 



the Kohn-Sham, Berry-phase, and Lagrange-multiplier 



-Dkk' = Im In det S^k' 



(24) 



and gi is the reciprocal lattice mesh vector in lattice di- 
rection i. (That is, k and k + g^ are neighboring k-points 

(i) 

in one of the iVj_ strings of k-points running in the recip- 
rocal lattice direction conjugate to a,.) Recall that the 
matrix of Bloch overlaps was defined in Eq. J3J. 

We now expand all quantities in orders of the perturba- 
tion, e.g., A(A) = A(°) + AAW+A 2 A( 2 ) + ..., etc. Similarly, 



we expand S kk ,(A) = 4k< + + \ 2 S& 

also define 



?(!) 



(2) 



k'k 



., and we 



(25) 



to be the inverse of the zero-order S matrix. Applying 
the 2n + 1 theorem to Eq. (|20(l , the variational second- 
order derivative of the total-energy functional is 



p(2) _ p(2) , p 



(2) 



F 



(2) 



LM 



(26) 



where 



F (2) _ 
^KS ~~ 



-. occ 

^EM> (0) 



"WlV-Vik. 



k.m 



+(V> 



mk I ^ext 



41 M 



(27) 



F (2) _ 
r BP — 



p(2) _ 
LM — 



ef ^ £ ■ at ^ (2 ) 
47r Z^ „«) Z^^k,k +gl 



(28) 



^-E [AgL«^ffi> + (^S» 



k,mn 



+A W 0) (V' (1 llV' ( l ) ) 



(29) 



In the Berry-phase term, Eq. I)28|l , we use the approach of 
Ref. flri! to obtain the expansion of lndetSkk' with respect 
to the perturbation. It then follows that 



D 



(2) 
kk' 



TmTr[2Sj£,Q 



; k'k 



k'k^kLWk'kl (30) 
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where , and Q are regarded as L x L matrices (L 
is the number of occupied bands), matrix products are 
implied, and Tr is a matrix trace running over the oc- 
cupied bands. Finally, in the Lagrange- multiplier term, 



Eq. (J2U, a contribution A 



(2) 

k,mn 



5 mn ) has 



been dropped from Eq. lj2"9"|) because the zero-order wave- 
functions, which have been calculated in advance, always 

satisfy the orthonormality constraints (V'mklV'ik) = &mn- 
Moreover, the zero-order Lagrange multipliers are made 
diagonal by a rotation among zero-order wavefunctions at 
each k point, and the first-order wavefunctions are made 
orthogonal to the zero-order ones at each iterative step, 
so that Eq. I|29|> simplifies further to become just 



Here, we have restored the notation e m k = A 
diagonal zero-order Lagrange multipliers. 



(0) 



(31) 



for the 



along the conjugate-gradient direction is easily de- 
termined to be 



1 dF& 

2 d6 



1=0 



d 2 F (2) 



d6 2 



}=0 



(36) 



3. Construction of the force- constant matrix 

To calculate phonon frequencies, we have to construct 
the force-constant matrix 



d 2 F 



dR va dR 



(37) 



Each diagonal element Q^p^p has already been obtained 
by minimizing the F^ in Eq. 126[) for the corresponding 
perturbation [ifi. The off-diagonal elements ^ va ^p can 
also be determined using only the first-order wavefunc- 
tions g using the (non-variational) expression 



2. Conjugate-gradient minimization 

The second-order expansion of the electric enthalpy 
functional in Eq. (|26[l is minimized with respect to 
the first-order wavefunctions using a "band-by-band" 
conjugate-gradient algorithmAii For a given point k and 
band m, the steepest-descent direction at iteration j is 
ICmkj) = dF^/d(u ( ^l I, where F^ is given by Eqs. $Z7\ 

Hp and (22). The derivatives of F$ and f[^ are 
straightforward; the new element in the presence of an 
electric field is the term 



dE, 



(2) 
BP 



ief 



£ ■ su 



'mk I 



(0 

i=l J, _L 



47r ~1 N 



|£>mk,k+ gi ) — |Ank,k- gi }) 

(32) 



where 



2W = ( \v$)Q Vk \v£?)Q Vk s£>,Q v1l ) 

\ / 1 



(33) 



In this equation, \u^)) and \u^,') are regarded as vectors 

of length L (e.g., Umk')' m = and vector-matrix leading to 

and matrix-matrix products of dimension L are implied 
inside the parentheses. Also, 



(0)\ 



c(i) _ /„(0) v , /..(l) I (0) v 

°kk',mri — \ u mkl"nk'/ ^ \ u mkl u nk'/ 



(34) 



are the first-order perturbed overlap matrices at neigh- 
boring k-points. The standard procedure translates the 
steepest-descent directions \Cmk.j) into preconditioned 
conjugate-gradient search directions \(p m k,j}- An im- 
proved wavefunction for iteration j + 1 is then obtained 
by letting 



^Hxc,fo;l u 'mk, ( ti l a 



20 f ( I (°) I (!) i (!) I (! 

+ <^Ll^a, W9 l«^>) dk + 5^.^/3 (38) 



where ua = dv cxt /dR„ a etc. 



B. Nonzero wavevector case 

In the case of a phonon of arbitrary wavevector q, the 
displacements of the atoms are essentially of the form 
bnva = A exp(iq ■ t„), where A is a complex number. 
However, a perturbation of this form does not lead by 
itself to a Hermitian perturbation of the Hamiltonian. 
This is unacceptable, because we want the second-order 
energy to remain real, so that it can be straightforwardly 
minimized. Thus, we follow the approach of Ref. and 
take the displacements to be 



bnua = Ae^"+Ae H 



(39) 



fext(r) = 4°t( r ) + AvexL Q , q (r) + A* v£> t Jr) 



+ aV 2) 



( n.q.q('*) A ^cxt,z^a,~q, — q(^") 
+ AA* Uext,i/a,q,-q( r ) + ^* A ^xt^Q,-q,q( r ) 



where 



(i) 

ext,i/a,±q 



d^ext(r) ±jq. t „ 



(40) 
(41) 



(1) X 
mk,j + l/ 



W mL') + %mk,i) 



(35) 



where 9 is a real number to be determined. Since the 
^-dependence of F^ 2 \9) is quadratic, the minimum of 



i-', 



,(2) 

ext,fa,±q,±q 



d 2 i; cxt (r) 



_ g ±iq.t» e ±iq-t n 



(42) 
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etc. Similarly, the field-dependent Bloch wavefunctions 
■0 and enthalpy functional F can also be expanded in 
terms of A and its hermitian conjugate as 

</W(r) = Vi°L(r) + A ^, q (r) + A* £ k) _ q (r) + ... (43) 
and 

F[£] = XF^[£]+F^[S\ + X*F^[£] 
+ A 2 F( 2 q [£] + 2AA*F( 2 l q [£] 



A* 2 F^_ q [£] 



(44) 



The first-order wavefunctions in response to a pertur- 
bation with wavevector q have translational properties 



Second, we consider the Berry-phase coupling term. 
The treatment of this term is rather subtle because, as 
mentioned above, the perturbed wavefunctions are now 
admixtures of parts with periodicity as in Eq. I|45l) and 
as in Eq. I|46l) , so that the usual Berry- phase formula for 
the polarization 9 cannot be used. A different approach 
is needed now in order to express the polarization in 
terms of the perturbed wavefunctions. For this purpose, 
we consider a virtual supercell in which the wavevectors 
k and q would be commensurate, and make use of the 
definition introduced by Resta^ specialized to the non- 
interacting case. The details of this treatment are de- 
ferred to the Appendix, but the results can be written in 
the relatively simple form 



^, q (r + R) = e < ( k +^- r ^, q (r) (45) 

that differ from those of the zero-order wavefunctions 

do) 



C(r + R)=e^l°»(r) 



(46) 



As a result, we cannot simply work in terms of perturbed 
Bloch functions or use the usual Berry-phase expression 
in terms of strings of Bloch functions. Also, in con- 
trast to the q=0 case, in which only one set of first-order 
wavefunctions was needed, we now need to solve for two 
sets V'mk ± q corresponding to the non-Hermitian pertur- 
bation at wavevector q and its Hermitian conjugate at 
wavevector q. 5 

We now proceed to write out the second-order en- 
ergy functional F^ [ip^; ?/£[ k ± ; £], corresponding to 
the sum of the quadratic terms in Eq. I|44(l , and minimize 
it simultaneously with respect to and . 

First, making the same decomposition as in Eq. (|26[1 . 
we find that the Kohn-Sham part is 



p(2) _ ef £ -Oj (2 ) 

i=l Jv l_ k 



(49) 



where 



-^k K&J — ±L [ a k,k+g^k+g,k J k,k+g-q 



(50) 



with Qk'k given by Eq. (|25ll and the superscript notation 
S {s,t) = d s+t S/d[x*) s dX t . From Eqs. Q and g5J}, we 
can write these explicitly as 



(1,0) 



°kk' 



(^°kle- 4g ' r |^, q ) 



(i) 



+ (^L_ q |e- ig - r |V^) . (51) 



(0) 



5, 



(0,1) 
kk'.r 



= (V> 



(0) 
k 



mk.q I ' 



l r nk , — q/ 



(52) 



P (2) _ p(2) r,/,(0) . ,/,(!) 1, P (2) r,/,(0).,(l) i ( A7 ~, 
^KS - ^q-qi^mk- V mk ,_ q J + ^-q.q[V mk , V mk , q J , (47) 

where 



0(1,1) 

kk',mn 



-igr 



+ (^ k ,-q|e- lg -|^,-q) • (53) 



(i) 



E 



(2) 



2n 



(27T) 



occ 

, /.._E 



BZ 



■ni 



b (0) \u {1) ) 

*mk,ql ^cxt.k+q^k+ql u mk,q/ 



(1) b (0) \u {1) ) 

mk.ql ^Hxc^k+q.k+ql^mk^q/ 
(1) L,(l) , ..(1) 



mk,ql cxt,k+q,k 



,(0) 



%xc,k+q,kl M mk/ 
(1) \ 



+ (u (0) \v il) +v {1) 

+ \%krcxt,k,k+q T ^Hxc.k.k+ql^mk.q/ 
+ (u (0) \v (2) \u (0) ))dk + -E {2) 

+ \ M mkl U ext,k,kl U mk/ I aK + ^ Hxc ' 



(48) 



(2) (2) 

Note that terms Eq t q and -E_ q! _ q vanish, essentially be- 
cause such terms transform like perturbations of wavevec- 
tor ±2q which, except when 2q equals a reciprocal lattice 
vector, are inconsistent with crystal periodicity and thus 
cannot appear in the energy expression. (If 2q is equal to 

(2) (2) 

a reciprocal lattice vector, _E q! q and still vanish, 

as can be shown using time-reversal symmetry.) 



Third, the treatment of the Lagrange-multiplier term 
is straightforward; in analogy with Eq. , we obtain 



F 



(2) 



LM 



-emk((^£k,ql 



^tq>+(V4i-qlV4t-q>) ■ (54) 



If we look closely at Eq. (|50|) . we see that the second 
term involves not simply pairs of k-points separated by 
the mesh vector g, but quartets of k-points, as illustrated 
in Fig. ^ Reading from left to right in the second term 
of Eq. H5Ufl , the k-point labels are k, then k + g — q, then 
k q, then k + g, and finally back to k. This is the loop 
illustrated in Fig.^ Each dark arrow represents a matrix 
element of S^'°\ S^, or Q; the gray arrow indicates 
the phonon q- vector. These loops arise because there are 
two kinds of coupling between k-points entering into the 
present theory. First, even in the absence of the phonon 
perturbation, wavevectors at neighboring k-points sepa- 
rated by mesh vector g are coupled by the £ ■ P term in 
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(a) 









| 




f 


< 

k 


i 



k+g-q 




F (2) 

remains in a quadratic form, the minimum of F<-V 
is again easily searched along the conjugate-gradient di- 
rection. Wavefunctions are updated over k-points one 
after another, and the first-order wavefunctions are up- 
dated. This procedure continues until the self-consistent 
potential is converged. Once the first-order responses of 
wavefunctions are obtained, the diagonal elements of the 
dynamical matrix are obtained by evaluating F^ 2 \ and 
the off-diagonal elements are obtained from a straight- 
forward generalization of Eq. I|38|) , 



FIG. 1: Pattern of couplings between k-points arising in (a) 
the first term, and (b) the second term, of Eq. I5U1 . Recipro- 
cal vector q is the phonon wavevector, while g is a primitive 
vector of the k-point mesh (indicated by thin horizontal and 
vertical lines). 



the energy functional. Second, the phonon introduces a 
perturbation at wavevector q. It is the interplay between 
these two types of inter-k-point coupling that is respon- 
sible for the appearance of these four-point loops in the 
expression for F^p . 

The implementation of the conjugate-gradient mini- 
mization algorithm proceeds in a manner very similar to 
that outlined in Sec. IIII A 21 Naively, one would have to 
work simultaneously with the two search-direction vec- 
tors 



|Cn*, q ) = dF^/diuHj 
iCmk.-q) = dF^/d(u^- . 



(55) 



where ± are the periodic parts of V'mk ±q- However, 
minimizing the second-order energy F^ with respect to 
two sets of first-order wavefunctions w„k,±q would double 
the computational cost and would involve substantial re- 
structuring of existing computer codes. We can avoid this 
by using the fact that the second-order energy is invari- 
ant under time reversal to eliminate one set of first-order 



wavefunctions ip^ _ q in favor of the other set V^k q 
lowing the approach given in Ref. [3. Specifically, the two 
sets of first-order wavefunctions are related by 



fol- 



</>i°kV) 
V^.qr) 



■Aw , 



r n _k, 



(56) 
(57) 



where #„k is an arbitrary phase independent of r. The 
arbitrary phase #„k cancels out in the expression of F^> 
since every term in F^ is independent of the phase of the 
first-order wavefunctions. Thus, we choose 6 n k = 1 for 
simplicity and write the second-order energy functional 
in terms of wavefunctions VVik,q only. 

The minimization procedure now proceeds in a man- 
ner similar to the zero-wavevector case, except that the 
calculation of the Berry-phase part involves some vector- 
matrix- matrix products as in Eq. (|33fl . but circulating 
around three of the sides of the loop in Fig. 2] Since 



2Q 



(2tt)3 



BZ 
(0) |„,(1) 



E 



(u {0) \v {1) 



(1) 



mk I ext,fa,k,k+q I mk,/x/3,q 



.(!) 



+ ( U mkl^Hxc, l 'a,k,k+ql M mk,^/3,q) 

+ ( u Lkl u oxt, t /a, W 3l M Lk) j^k 



(2) 

Hxc,^a,/i./3 " 



(58) 



IV. 



TEST CALCULATIONS FOR III-V 
SEMICONDUCTORS 



In order to test our method, we have carried out cal- 
culations of the frequency shifts induced by electric fields 
in two III-V semiconductors, AlAs and GaAs. We have 
chosen these two materials because they are well-studied 
systems both experimentally and theoretically, and be- 
cause the symmetry allows some phonon mode frequen- 
cies to shift linearly with electric field while others shift 
quadratically. Since our main purpose is to check the in- 
ternal consistency of our theoretical approach, we focus 
on making comparisons between the shifts calculated us- 
ing our new linear-response method and those calculated 
using standard finite-difference methods. Moreover, as 
mentioned at the start of Sec. Ill B II we have chosen to 
neglect changes in phonon frequencies that enter through 
the electric-field induced strains (piezoelectric and elec- 
trostrictive effects), and we do this consistently in both 
the linear-response and finite-difference calculations. For 
this reason, our results are not immediately suitable for 
comparison with experimental measurements. 

Our calculations are carried out using a plane-wave 
pseudopotential approach to density-functional theory. 
We use the ABINIT code package^ which incorpo- 
rates the finite electric field method of Souza et al£ for 
the ground-state and frozen-phonon calculations in finite 
electric field. We then carried out the linear-response cal- 
culations with a version of the code that we have modified 
to implement the linear-response formulas of the previous 
section. 

The details of the calculations are as follows. We use 
Troullier-Martins norm-conserving pseudopotentials, 20 
the Teter Pade parameterization^! of the local-density 
approximation, and a plane- wave cutoff of 16 Hartree. A 
10x10x10 Monkhorst-Pack 2 ^ k-point sampling was used, 
and we chose lattice constants of 10.62 A and 10.30 A for 
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TABLE I: Calculated frequency shifts, in cm" 1 , induced by 
an electric field of 5.14 x 10 8 V/m applied along x in GaAs 
and AlAs. 'FD' are the results of finite-difference (frozen- 
phonon) calculations in which atoms are displaced by hand 
and restoring forces are calculated, while 'LR' refers to the use 
of the linear-response developed here. The L and X points are 
at (2ir/a)(l, 1, 1) and (27r/a)(l, 0, 0) respectively. 





GaAs 




AlAs 




Mode 


FD 


LR 


FD 


LR 


r oi a 


-3.856 


-3.856 


-5.941 


-5.941 


r 02 a 


-0.282 


-0.281 


-0.300 


-0.299 


r 03 a 


3.548 


3.548 


5.647 


5.647 


L LO 


2.701 


2.703 


4.282 


4.282 


L TOl 


-3.749 


-3.749 


-5.663 


-5.663 


L T02 


0.567 


0.564 


0.952 


0.952 


X LO 


0.050 


0.050 


-0.243 


-0.243 


X TOl 


-3.953 


-3.953 


-6.083 


-6.083 


X T02 


3.753 


3.753 


5.919 


5.919 



"The non-analytic long-range Coulomb contributions are excluded 
for the r modes. 



AlAs and GaAs, respectively. The crystals are oriented 
so that the vector (a/2)(l, 1, 1) points from a Ga or Al 
atom to an As atom. 

Table [I] shows the changes in phonon frequencies re- 
sulting from an electric field applied along a Cartesian 
direction at several high-symmetry q-points in GaAs 
and AlAs. Both the electronic and ionic contributions, 
Eqs. (|7I8|I . are included. We first relaxed the atomic co- 
ordinates in the finite electric field until the maximum 
force on any atom was less than 1CP 6 Hartree/Bohr. We 
then carried out the linear-response calculation, and in 
addition, to check the internal consistency of our linear- 
response method, we carried out a corresponding calcu- 
lation using a finite-difference frozen-phonon approach. 
For the latter, the atoms were displaced according to 
the normal modes obtained from our linear-response 
calculation, with the largest displacement being 0.0025 
Bohr. (Because the electric field lowers the symmetry, 
the symmetry-reduced set of k-points is not the same 
as in the absence of the electric field.) The agreement 
between the finite-different approach and the new linear- 
response implementation can be seen to be excellent, with 
the small differences visible for some modes being at- 
tributable to truncation in the finite-difference formula 
and the finite density of the k-point mesh. 

In Table UTI we decompose the frequency shifts into the 
ionic contribution A0; on (q; £ ) and the electronic contri- 
bution A0 G i(q;£) defined by Eqs. © and (7J, respec- 
tively, calculated using the linear-response approach. It 
is clear that the largest contributions are ionic in origin. 
For example, the large, roughly equal and opposite shifts 
of the 01 and 03 modes at T arise from the ionic terms. 
However, there are special cases (e.g., 02 at T and LO at 
X) for which the ionic contribution happens to be small, 
so that the electronic contribution is comparable in mag- 
nitude. 



TABLE II: Same as in Table but with the frequency shifts 
decomposed into ionic and electronic contributions as defined 
in Eqs. (JHJ and J7J respectively. 





GaAs 




AlAs 






Ion 


Elec. 


Ion 


Elec. 


r oi a 


—3.659 


—0.198 


—5.684 


—0.257 


r 02 a 


-0.146 


-0.135 


-0.123 


-0.177 


r 03 a 


3.655 


-0.107 


5.589 


0.058 


L LO 


2.341 


0.362 


3.633 


0.649 


L TOl 


-3.486 


-0.262 


-5.628 


-0.034 


L TO2 


1.181 


-0.617 


1.658 


-0.707 


X LO 


0.122 


-0.073 


-0.033 


-0.209 


X TOl 


-3.411 


-0.543 


-5.658 


-0.424 


X T02 


3.388 


0.365 


5.609 


0.310 



"The non-analytic long-range Coulomb contributions are excluded 
for the r modes. 



The pattern of ionic splittings appearing at T can be 
understood as follows. Because the non-analytic long- 
range Coulomb contribution is not included, the three 
optical modes at T are initially degenerate with frequency 
luq in the unperturbed lattice. A first-order electric field 
along x induces a first-order relative displacement u x of 
the two sublattices, also along x. By symmetry consid- 
erations, the perturbed dynamical matrix is given, up to 
quadratic order in u x , as 

( l + Hul \ 

D(T) = ui%\ 1 + vul ku x . (59) 

\ KU X 1 + VV? X J 

The off-diagonal n term arises from the E xyz coupling 
in the expansion of the total energy in displacements; 
this is the only third-order term allowed by symmetry. 
The \i and v terms arise from fourth-order couplings 
of the form E xxxx and E xxyy respectively. The eigen- 
values of this matrix are proportional to 1 + fiu^ and 
1 ± ku x + vu x . Thus, two of the modes should be per- 
turbed at first order in the field-induced displacements 
with a pattern of equal and opposite frequency shifts, 
while all three modes should have smaller shifts arising 
from the quadratic terms. This is just what is observed in 
the pattern of frequency shifts shown in Table [n] (The 
symmetry of the pattern of electronic splittings is the 
same, but it turns out that the linear shift is much smaller 
in this case, so that for the chosen electric field, the linear 
and quadratic contributions to the electronic frequency 
shift have similar magnitudes.) A similar analysis can 
be used to understand the patterns of frequency shifts at 
the L and X points. 

We have also plotted, in Fig. [3 the calculated total 
frequency shift A0 o i(q) + A0j on (q) and its electronic con- 
tribution A0 o i(q) along the line from T to L for the case 
of AlAs. (The 'LO' and 'TO' symmetry labels are not 
strictly appropriate here because the electric field along 
x mixes the mode eigenvectors; the notation indicates 
the mode that would be arrived at by turning off the 
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FIG. 2: Frequency shifts induced by an electric field of 
5.14xl0 8 V/m along x in AlAs, plotted along F to L. Filled 
and open symbols indicate the total shift A0 c i + A0i on and 
the electronic contribution A0 c i respectively. 



FIG. 3: Frequencies of LO and TO modes at L in AlAs as a 
function of electric field (where 10~ 3 a.u. = 5.14 x 10 8 V/m) 
applied along x. The symbols have the same interpretation 
as in Fig. |2] 



field.) In contrast to the results presented in Tables ITITT1 
the frequencies at T in Fig. [3 were computed by includ- 
ing the long-range non-analytic Coulomb contribution for 
q || (111) in order to extend the curves to q = 0. (Because 
the direct linear-response calculation of the dynamical ef- 
fective charge and dielectric susceptibility tensors are not 
yet developed and implemented in the presence of a finite 
electric field, the needed tensor elements were computed 
by finite differences.) It is clearly evident that the elec- 
tronic terms remain much smaller than the ionic ones for 
all three optical modes over the entire branch in (/-space. 

Returning now to the comparison between our ex- 
act theory of Sec. Ill B II and the approximate theory of 
Sec. Ill B 21 we compare the equilibrium positions and 
phonon frequencies predicted by these theories in Ta- 
ble IIIII Recall that Rs is calculated in the approxi- 
mate theory by using Eq. (|10fl . Using this force, the 
ion coordinates were again relaxed to a tolerance of 10 -6 
(Hartree/Bohr) on the forces. It can be seen that Rg 
is predicted quite well by the approximate theory, with 
errors of only ~2%, confirming that the displacements 



TABLE III: Comparison of ionic displacements and frequency 
shifts at the L point in GaAs as computed by the approximate 
and exact approaches of Sec. Ill B 21 and IHB II respectively, 
again for an electric field of 5.14 x 10 8 V/m along x. Re is the 
induced displacement of the cation sublattice along x, and 
the A0i on are ionic contributions to the frequency shifts as 
defined in Eq. @. 







Re 


A0 


on(L) (cm" 


X ) 






(10" 3 A) 


LO 


TOl 


T02 


GaAs 


Approx. 


5.07 


2.63 


-3.89 


1.37 




Exact 


4.95 


2.34 


-3.49 


1.18 


AlAs 


Approx. 


5.69 


3.75 


-5.66 


1.65 




Exact 


5.62 


3.63 


-5.63 


1.66 



can be calculated to good accuracy using a linearized 
theory for this magnitude of electric field. The changes 
in the phonon frequencies resulting from these displace- 
ments (evaluated at zero and non-zero field for the ap- 
proximate and exact theories respectively) are listed in 
the remaining columns of Tabic IIIII The discrepancies 
in the phonon frequencies are now somewhat larger, ap- 
proaching 15% in some cases. This indicates that the 
approximate theory is able to give a moderately good 
description of the phonon frequency shifts of GaAs in 
this field range, but the exact theory is needed for ac- 
curate predictions. (Also, recall that the approximate 
theory does not provide any estimate for the electronic 
contributions, which are not included in Table ITTTl 'l 

Finally, we illustrate our ability to calculate the non- 
linear field dependence of the phonon frequencies by pre- 
senting the calculated optical L-point phonon frequencies 
of AlAs in Fig. |31 as a function of electric field along x. 
These are again the results of our exact theory, obtained 
by including both ionic and electronic contributions. The 
two TO modes are degenerate at zero field, as they should 
be. All three modes show a linear component that dom- 
inates their behavior in this range of fields. However, 
a quadratic component is also clearly evident, illustrat- 
ing the ability of the present approach to describe such 
nonlinear behavior. 



V. SUMMARY AND DISCUSSION 

We have developed a method for computing the 
phonon frequencies of an insulator in the presence of 
a homogeneous, static electric field. The extension of 
density-functional perturbation theory to this case has 
been accomplished by carrying out a careful expansion 
of the field-dependent energy functional £ks + Q£ ■ P, 
where P is the Berry-phase polarization, with respect to 
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phonon modes both at q = and at arbitrary q. In 
the general case of nonzero q, there is a subtle interplay 
between the couplings between neighboring k-points in- 
troduced by the electric field and the further-neighbor 
couplings introduced by the q- vector, so that terms arise 
that require the evaluation of four-sided loops in k-space. 
However, with the judicious use of time-reversal symme- 
try, the needed evaluations can be reduced to a form that 
is not difficult to implement in an existing DFPT code. 

We have carried out test calculations on two III-V 
semiconductors, AlAs and GaAs, in order to test the cor- 
rectness of our implementation. A comparison of the re- 
sults of linear-response and finite-difference calculations 
shows excellent agreement, thus validating our approach. 
We also decompose the frequency shifts into "lattice" 
and "electronic" contributions and quantify these, and 
we find that the lattice contributions (i.e., those resulting 
from induced displacements in the reference equilibrium 
structure) are usually, but not always, dominant. We 
also evaluated the accuracy of an approximate method for 
computing the lattice contribution, in which only zero- 
field inputs are needed. We found that this approximate 
approach gives a good rough description, but that the 
full method is needed for an accurate calculation. 

Our linear-response method has the same advantages, 
relative to the finite-difference approach, as in zero elec- 
tric field. Even for a phonon at T, our approach is more 
direct and simplifies the calculation of the phonon fre- 
quencies. However, its real advantage is realized for 
phonons at arbitrary q, because the frequency can still be 
obtained efficiently from a calculation on a single unit cell 
without the need for imposing commensurability of the q- 
vector and computing the mode frequencies for the corre- 
sponding supercell. We also emphasize that the method 
is not limited to infinitesimal electric fields. We thus ex- 
pect the method will prove broadly useful for the study of 
linear and nonlinear effects of electric bias on the lattice 
vibrational properties of insulating materials. 
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APPENDIX 

The formula for the electric polarization given in the 
original work of King-Smith and Vanderbilt 9 is not a suit- 
able starting point for the phonon perturbation analysis 
that we wish to derive here, because a perturbation of 
nonzero wavevector q acting on a Bloch function gener- 
ates a wavefunction that is no longer of Bloch form. That 
is, while the zero-order wavefunction V'mk transforms as 
e ik R urK j er a translation by R, the first-order wavefunc- 
tion tp^. q transforms as e l ( k+q )' R . 



To solve this problem, we first restrict ourselves to the 
case of a regular mesh of N = N\ X JVjj X N3 k-points in 
the Brillouin zone. As is well known, one can regard the 
Bloch functions at these k-points as being the solutions 
at a single k-point of the downfolded Brillouin zone of an 
N% x N2 x N3 supercell. Then, as long as the wavevec- 
tor q is a reciprocal lattice vector of the supercell, or 
q = mibx/iVi + rri2i>2/N2 + m-sh^/Ns, the phonon per- 
turbation will be commensurate with the supercell, and 
the perturbed wavefunction will continue to be a zone- 
center Bloch function of the supercell. We thus restrict 
our analysis to this case. 

A formula for the Berry-phase polarization for single- 
k-point sampling of a supercell has been provided by 
Resta^ Starting from a general many-body formulation 
in terms of a definition of the position operator suitable 
for periodic boundary conditions, and then specializing to 
the case of a single-particle Hamiltonian, Rcsta's deriva- 
tion leads to 

P = ^E7.a a (60) 

where the Berry phase in lattice direction a is given by 
7« = -Im In det M a . (61) 

Here 

M a , ss , = (^ s |e- ig - r |Vv> (62) 

where g a = h a /N a is the primitive reciprocal mesh vec- 
tor in lattice direction a and s runs over all of the occu- 
pied states of the supercell. Expanding the matrix M a 
in powers of A and A*, 

M a (X, A*) = M <°>°) + \m£>V + \*Mi°>V + A 2 M( 2 >°) 
+ \X\ 2 M^ +X* 2 M^ + ... , (63) 

the expansion of In det M a (X, A*) takes the formic 

lndetM a = lndetM^°' 0) 

+ XTv[M^Q a ] + X* Tr[Mi°VQ a ] 
+ A 2 Tr[2A/( 2 <°>Q Q - M^Q a M^Q a ] 
+ A* 2 Tr[2M(°' 2 )g a - Mt 1] Q a M^Q a ] 
+ \X\ 2 Tr[2M^Q a - M^Q a M^Q a 

-M^Q a M^Q a ] 
+ higer order terms , (64) 

where Q a = [M&°' 0) ]- 1 . 

From the physical point of view, the terms proportional 
to A, A*, A 2 , and A* 2 should vanish as a result of transla- 
tional symmetry. For example, a term linear in A should 
transform like e lq R under translation by a lattice vector 
R, but such a form is inappropriate in an expression for 
the energy, which must be an invariant under transla- 
tion. We have confirmed this by explicitly carrying out 
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the matrix multiplications for these terms and checking 
that the traces are zero. Using the cyclic property of the 
trace to combine the last two terms, we find that the 
overall second-order change in In det M a is 

(lndetA/ Q )( 2 ) = 2|A| 2 Tr[M< [ 1 < 1 )Q Q 

-Mt 0) Q a M^Q a ] . (65) 

In our case, the orbitals ip s appearing in Eq. 1|62[1 are 
the perturbed wavefunctions originating from the unper- 
turbed states labeled by band m and k-point k of the 
primitive cell, so that we can let tf> s — > ipmk and 

M a , m k,m>k> = (V'mk|e~ 4Sa ' r |VWk'} • (66) 
Substituting Eq. into Eq. we find 
r(i,o) 



M 



Q,mk,m'k' 



M (o,i) 

a.rnk.m'k' 



(0)l — id*> 



mk I 



(i) 

mk,— q 
(0) 



1&> . (67) 



mk.q I ' 



Id,-,) 

k' 



(68) 



M 



Cbi) 

a.mk.m'k' 



(4> {1 l |e- ,;g °- r |V' (1 L ) 

\ T^mk.a 1 1 ^m'k .a/ 



mk,q 

/,(! 



+ (VW,-ql 

and Q Q = [Af! 0,0 ^] -1 where 



m'k y ,q/ 

Id,-,), (69) 



71/ 



(0,0) 



= ^ (0) le- is °- r l^ (0) ) 

o.mk.m'k' \ ^mk I I r m'k'/ * 



(70) 



The transformation properties of the zero- and 
first-order wavefunctions under translations, given by 
Eqs. l (4*S|l and jlSjl. impose sharp constraints upon which 
of the terms in Eqs. I|67I70|I can be non-zero. For exam- 
pie, for Afi 1,0) in Eq. (£3, the term (V&^'IV^k'.q) 
is only non-zero if k = k' + q — g Q . Similarly, Q a ,mk.m'k' 
is only non-zero if k = k' + g Q . In practice, we define 
primitive-cell-periodic functions 



(71) 



and 



so that 



^L, q (r) = e- <(k+q) - r ^l q (r) 



(72) 



M 



(i,o) 

a,mk,m'k' 



(1,0) 



cl- 1 -, 1 



kk'.mm' °k,k'+q-g Q 



M (0,l) _ o(0,l) r 

J ' J a,mk,ro'k' — D kk', mm' °k,k'-q- gcv , 

)if(l,l) _ c(l,l) I 

M a,mk,m'k' — D kk'. mm' °k,k'-g Q 7 



Jw a,mk,m'k' ~ °kk\mm' u k,k'-g Q 



-,(0,0) 



where 



(73) 



c(i-o) 

kk' , mm' 


\ U mk 


U m'k' 


.q/ ' \"mk,-q 


|t&> 


9(0,1) 

kk' , mm' 


\ U mk 


U m'k' 


) + (u (1) 

.-q/ T \"mk,q 


It&> 


0(1,1) 

kk' , mm' 


\ U mk, 


ql m' 


k',q/ T \"mk,- 


ql"m'k' 


0(0,0) 

kk',mm' 


( U (0) 
\"mk 


(0) 
U m'k' 







(74) 



(subscript a is now implicit). Defining Qk'k = [S^? ] _1 
and taking into account the constraints on k-points em- 
bodied in the delta functions in Eq. 1)73(1 . the two terms 
in Eq. ljBl)j) become 



Tr[M£>VQ c 



E Tr ^k!kigQk + g,k] (75) 



and 



Tr[M^Q a M^Q a ] = £ Tr[s£ k °j g _ q 



X Qk+g-q.k-qS'k-q.k+gQk+g.k 



(76) 



In these equations, the trace on the left-hand side is over 
all occupied states of the supercell, while on the right- 
hand side it is over bands of the primitive cell. These are 
the terms that appear in Eq. (|50|l in the main text, and 
that determine the pattern of k-point loops illustrated in 
Fig. ID 
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